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Two weak coupling Continuous Time Quantum Monte Carlo (CTQMC) methods are shown to 
be equivalent for Hubbard-like interactions. A relation between these CTQMC methods and the 
Hirsch-Fye Quantum Monte Carlo (HFQMC) method is established, identifying the latter as an 
approximation within CTQMC and providing a diagrammatic interpretation. Both HFQMC and 
CTQMC are shown to be equivalent when the number of time slices in HFQMC becomes infinite, 
implying the same degree of fermion sign problem in this limit. 
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PACS numbers: 02.70.Ss,71.10.Fd,71.27.+a 

Introduction. Hirsch-Fye Quantum Monte Carlo 
(HFQMC) is a standard method for the simulation of 
quantum lattice models [3, 0, B 0, 0|- However, dur- 
ing the past decade, new QMC methods have emerged, 
which are based on stochastic sampling of diagrams in 
a perturbative expansion These new methods avoid 
systematic errors due to finite discretization in the imag- 
inary time, and are commonly referred to as "Continuous 
Time QMC" (CTQMC). With new variants of CTQMC 
appearing, a comparison of these formalisms becomes 
important. In this work, we consider two seemingly 
different CTQMC methods proposed by Rombouts [7[ 
and Rubtsov Q based on the expansion of the inter- 
action term of Hamiltonian in the perturbation series, 
also known as "weak coupling" CTQMC. We show that 
for Hubbard-like interactions these methods are equiva- 
lent. We also show that HFQMC can be interpreted as 
a summation of a specific subset of diagrams present in 
CTQMC. 

The equivalence of CTQMC algorithms by Rombouts 
and Rubtsov. These methods consider the perturbative 
expansion of the partition function in powers of the in- 
teraction and then sample the resulting series of multi- 
dimensional integrals stochastically. We will use a path 
integral formalism to illustrate this. Here, the partition 
function is written as an integral over the Grassman vari- 
ables 17,7?*: Z = J Vn*Vne' s{r1 ' , with the action 



S(r)*, V ) = So(r,*,r,)- dTV(r,*(r), V (r)) , (I) 
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sion we consider Hubbard like repulsive interaction 
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In the Rombouts method 0], a constant K is intro- 
duced to shift the reference free energy and the resulting 
series expansion for the partition function can be written 
as: 
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The following identity is then used to decouple the inter- 
action terms and introduce an auxiliary field s: 
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where cosh 7 = 1 + . The resulting series for the 



partition function is: 



Z = e 



-K 



,-So 



E 

kjrs 



( K 



(5) 



x e 



where multiple sums and integrals are denoted as: 



where So is the bare part of S, and V is the interacting 
part of the Hamiltonian H. For the purposes of discus- 
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The fermion degrees of freedom can now be integrated 
out, and the partition function can be rewritten as [Io| : 
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where G\ is the Green's function for a particular con- 
figuration of auxiliary fields, and is related to the non- 
interacting Green's function Q® by a Dyson's equation: 



-fW, 
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with W} Sl> = diag(o-s i ) and [Q°] pq = G°(j p , t p ; j q , r q ) 
being k x k matrices. Finally, QMC is used to perform 
the multidimensional sum (Eq. [6]) over different expan- 
sion orders and configurations of the auxiliary fields. For 
this, a Markov process is set up that samples the configu- 
rations of random auxiliary fields {s^} with weight given 
by the product of determinants in Eq. [7] 

In the Rubtsov method 
rewritten as: 
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thereby introducing auxiliary fields s. This amounts to 
introducing a shift in the free energy: 



K = [3UN C 



(10) 



The auxiliary fields s suppress the oscillating sign of the 
integrand in the perturbative expansion Q: 
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The fermion degrees of freedom can be integrated out, 
and the partition function becomes [ill ]: 
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Again, the product of determinants in Eq. [12] gives the 
weight in QMC to evaluate the multidimensional sum 
over the configurations of auxiliary fields {§i}- 

We now show that the two expansions (O and (fT2|) are 
equivalent (term by term) and that the auxiliary fields 
{si} and {§i} are equivalent as well. Using Eq. the 



inverse Green's function G 
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Using this, and the fact that 
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from which we deduce that both algorithms are equiva- 
lent if a = a* , which is the same as requiring that Eq. [TU] 
holds for freely adjustable parameters K and a in these 
methods. Both algorithms must have the same degree 
of sign problem and statistics of measurements (such as 
auto-correlation time), as long as the above mentioned 
condition for the parameters K and a is satisfied. 

The relation between HFQMC and CTQMC. The 
derivation of the Hirsch-Fye algorithm involves breaking 
up the partition function using a Trotter decomposition 
and decoupling the quartic part of the Hamiltonian with 
the transformation [2L Il2|: 



,-A-r£/[ntni-3(nt-fni)] = I 
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where cosh A = e ArU ^ 2 . The resulting partition function 
takes the well-known form [l|, 0] : 
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where G^ 3 ^ is the Green's function for a particular con- 



(11) figuration of auxiliary fields {sj}: 
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The product of determinants (|17[) yields a sampling 
weight for a corresponding configuration of the auxiliary 
fields {sj}. It is very similar to the sampling weight in 
CTQMC (Eq. [7J, as are the transformations employed 
(Eqs.rR)1andHl) and the update formulas (Ea.[T51and[5|l. 
The only difference is that in HFQMC, the number of the 
auxiliary fields is fixed to /chf = (3N c /At, and they are 
distributed evenly in the imaginary time. In addition, 
the parameter A plays the same role as the parameter 7 
in CTQMC to couple the auxiliary fields to the fermion 
spin. In fact, one can formulate a set of restrictions, un- 
der which CTQMC reduces to HFQMC: 

1. Restrict the expansion order k in CTQMC equal 
to the number of the auxiliary fields /sjjf in the 
HFQMC and distribute them evenly in the imagi- 
nary time interval (0 . . . [3). 
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2. Set the strength of the auxiliary field in CTQMC: 
7 = A. In terms of CTQMC parameters K or a, 



this condition is equivalent to: 

0UN c 0UN c 
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3. Restrict the Monte-Carlo moves to flipping the aux- 
iliary fields associated with the interaction vertices; 
shifting vertices in imaginary time is not allowed. 

These restrictions imply that only a subset of diagrams 
with fixed expansion order and equidistant auxiliary 
fields are sampled in HFQMC, whereas in CTQMC, all 
diagrams of variable order and all possible sets of auxil- 
iary field configurations contribute (see Fig. [T]) . 



Since the attractive Hubbard model has no sign problem, 
the parameter a can be set equal to zero. However, for 
a > 0, the relation to HFQMC is again given by the same 
set of restrictions as defined above (including Eqs. 119120)1 . 
Similarly, the relation between CTQMC and HFQMC is 
preserved in case of nonlocal density-density interactions. 

Small At limit. When At — » 0, systematic errors 
in HFQMC arc eliminated and in this sense HFQMC 
and CTQMC are equivalent. The relationship described 
above will also hold for a — * oo (see Eq. [2U1) . In the 
discussion above, HFQMC is interpreted as sampling just 
one order in series expansion. To understand this, we 
need to revisit the sampling and measurement procedure 
in the CTQMC. The expectation value of any operator 
can be written as a series expansion: 
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FIG. 1: Example of the diagrams sampled in HFQMC (top) 
and CTQMC (bottom). The light and dark lines denote prop- 
agation of spin-up and spin-down fermions, respectively. The 
block arrows represent the auxiliary field, associated with an 
interaction vertex. These vertices are distributed evenly in 
HFQMC, with interval At, while the distribution, positions 
and the number of vertices are arbitrary in the CTQMC. Note 
that both connected and disconnected diagrams are sampled 
in both methods. In fact, product of determinants (Eqs. [121 
and !17[) accounts for all (fe!) 2 possible diagrams for a partic- 
ular set of vertices. 



The relation between CTQMC and HFQMC persists 
for the attractive Hubbard model (U < 0). Here, a dif- 
ferent discrete Hubbard Stratonovich transformation is 
used: 
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with cosh A = e Ar '' 7 '/ 2 . The corresponding form for the 
interaction in CTQMC also has to be modified: 
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In both variants of CTQMC the evaluation of this sum 
is done with importance sampling, and the weight (or 
the "guiding function") is taken to be equal to the corre- 
sponding contribution to the partition function Z^, with 
Qk — being the Monte Carlo estimator for a fixed 
order of expansion. Of course, Z^ depends on the con- 
figuration of the auxiliary fields {s^}, so the actual esti- 



mator is gj. 



'k 



However, for this discussion, we 



are only interested in how this estimator depends on the 
expansion order, so we assume that the auxiliary fields 
are already summed. 

The series expansion for the partition function (Eq. llip 
defines a distribution (see Fig. [2]) with mean value [11] : 



dr (V(t)) 



(24) 



This can be generalized for higher factorial moments: 



((k) n ) z ={k(k-l)...(k-n+l)) z 
= (-1)" fdn...J dT„ (T T V( n ) . . . V(t„)} 
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For the Hubbard model with sufficiently large a, these 
moments scale as: 



lim {(k) n ) z = ([3UN c a 2 ) n =p n 



(26) 



which is a property of Poisson distribution P p (k) = — 
with parameter p — f3UN c a 2 . Of course, for large p, the 
Poisson distribution approximates a normal distribution 
(see Fig. [2]). In a similar way as the series expansion for 
the partition function defines its distribution, the expan- 
sion for the Green's function (or any measurable quan- 
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FIG. 2: The distribution of the contributions Zk (Gk) to the 
partition function (Green's function) vs. the expansion order 
k for the 2D Hubbard model with N c = 16 sites, U = W 
(bandwidth), filling / = 0.75, f3 = 4.5, a = 1.01. Note 
that the distributions almost perfectly overlap, and their ra- 
tio (proportional to gk) varies very little with k. Restricting 
diagrams to one order k — &hf and restricting vertices to dis- 
crete imaginary time grid results in the Hirsh-Fye algorithm 
with Art = 0.125. 
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FIG. 3: Comparison of the average sign for the 2D Hubbard 
model with N c = 24E, U = St, filling= 0.95, f3t = 4.5. In 
CTQMC, the average sign depends on parameter a (lower 
horizontal axis), while in HFQMC, it depends on the size of 
the time slice Art (upper horizontal axis). Marks on both 
axes are related by Eq. [20] 



tity) defines another distribution: 

G{n,r 3 ) = J2G k (r l ,r J )=J2 { - Z ^ f dr x ... 

k k ' ^° 

r/3 

... dT k ■ {T T c{T i )c\T j )V{T 1 )...V{T k )) , (27) 
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which is characterized by its factorial moments: 

(-1)™ f p f 
(( fc )n>G„ = 777 7 / d n ... dr n x 

x (Trcin^iT^Vin) . . . V{r n )) . (28) 

In general, this distribution is different from the one de- 
fined by the expansion of the partition function. How- 
ever, in the limit when a — * oo, the factorial moments of 
Green's function distribution scale as: 

lim <(fc)„) G = (l3UN c a 2 ) n = lim ((k)„) z . (29) 

Since all the moments for both distributions are the same, 
the distributions are the same as well in this limit, and 
the estimator g k becomes a constant, independent of the 
expansion order k (see Fig. [5]). Thus, the sum over all 
expansion orders k can be replaced by any single term 
corresponding to a fixed value of k = /chf ■ That explains 
why sampling just one single order in the expansion for 
the partition function (as is done in HFQMC) gives the 
same exact result when At — > 0. 

Computational implications. When the product of de- 
terminants (Eqs. I7I12I17|) is not positive definite, its ab- 
solute value is taken as a weight in QMC. This approach 



fails, if the average sign of the product of determinants 
becomes small. This is the infamous fermion sign prob- 
lem, the main limitation in any fermion QMC method. 
From the discussion above, it follows that both HFQMC 
and CTQMC have the same degree of sign problem when 
a — > oo and At — > 0. For typical finite values of At and 
a, the difference in average sign is still small (see Fig. [3]) 
and depends on model parameters. Altogether, we find 
that neither of the methods has a definite advantage in 
terms of the degree of the sign problem. Also, the auxil- 
iary fields enter the same way in both methods, and corre- 
lations in these fields give information about the spin and 
charge correlations in the repulsive and attractive Hub- 
bard models, respectively. Thus, optimization strategies 
developed for HFQMC can be applied to CTQMC. 

Conclusions. We have investigated two weak coupling 
CTQMC methods proposed by Rombouts and Rubtsov, 
and shown that they are equivalent for a certain choice 
of freely adjustable parameters in these methods. We 
also established the relation between the CTQMC meth- 
ods and HFQMC method and identified the latter as an 
approximation within CTQMC where the Monte Carlo 
sum is restricted to a certain subset of diagrams. We 
have shown that this approximation becomes exact in 
the limit when an infinite number of time slices is taken 
in HFQMC, implying that both methods have the same 
degree of the sign problem in this limit. 
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